function [psi,phi] = cmsc3(psi,phi,Q,in,mu,w,dx,g)
% function [psi,phi] = cmsc3(psi,phi,Q,in,mu,w,dx)
%   Inputs:
%           psi     -- fine mesh angular flux
%           phi     -- fine mesh scalar flux
%           Q       -- fine mesh external source
%           in      -- inpute structure
%           mtt     -- material index for each fine mesh
%           mu      -- angle set
%           wt      -- angle weight
%           dx      -- fine mesh delta x
% 
%   This function is a different form of CMSC.  Instead of using 
%   discontinuity factors w/r to currents at coarse mesh edges, we
%   introduce a perturbed source term that begins a coarse iteration
%   with the right angular fluxes at the the edges.  Again, we cast
%   this in a S2-like framework.

fixed = 0; % fixed perturbation approach

% allocate
sigt  = zeros( length(in.xfm), 1 );   % homogenized total cross-section
sigs  = zeros( length(in.xfm), 1 );   % homogenized scatter cross-section
d_sigs  = zeros( length(in.xfm), 2 ); % perturbed scatter cross-section, one per coarse direction
s_c   = zeros( length(in.xfm), 1 );   % volume-averaged external source
phi_c = zeros( length(in.xfm), 1 );   % coarse mesh scalar flux
psi_c = zeros( length(in.xcm), 2 );   % coarse mesh edge angular flux
psi_c_ref = zeros( length(in.xcm), 2 );   % reference coarse mesh edge angular flux
Jp    = zeros( length(in.xcm), 1 );   % coarse mesh edge rightward (NOT J+ as in description)
Jm    = zeros( length(in.xcm), 1 );   % coarse mesh edge leftward (NOT J- as in description)
Jref  = zeros( length(in.xfm), 2 );   % reference coarse mesh total incoming and outgoing partial current
Jhom  = Jref;                         % homogenized coarse mesh total p cur
q_c   = zeros( length(in.xfm), 1 );   % coarse mesh emission density
Jp_c = Jp; 
Jm_c = Jm;
% coarse mesh delta x's
h = in.xcm(2:end)-in.xcm(1:end-1);

s = zeros(sum(in.xfm),1);  % the angle-integrated ext src (i.e ausume isotropic)
for i = 1:length(s)
    s(i) = sum( Q(i,:,1)'.*w(:) );
end

% The partial current is defined e.g. 
%    Jp = int_{mu>0} mu*psi dmu = sum w_i mu_i psi_i, for mu_i > 0
% What we'll do here is define the partial currents using whatever the fine
% mesh quadrature is (S2, S8, etc.), compute the partial currents, and then
% simply use the S2-like approach to define psi_c = Jp/mu_c.  We could
% stick in the full angular space, an option to flesh out later.

% using an S2-like quadrature
mu_c = 0.5773502691;
wt_c = 1;

Jp(1) = partcur( 1, psi, 1, mu, w, in );
Jm(1) = partcur(-1, psi, 1, mu, w, in );
for i = 1:length(in.xfm)  
    % note, the homogenization is meaningless until I implement a way to
    % specify coarse meshes that contain heterogeneities
    idx1     = 1 + sum( in.xfm(1:(i-1)) );            % lower index
    idx2     = sum( in.xfm(1:(i  )) );                % upper index   
    Jp(i+1)  = partcur( 1, psi, idx2+1, mu, w, in );  % J+ -->
    Jm(i+1)  = partcur(-1, psi, idx2+1, mu, w, in );  % J- <--
    phi_c(i) = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1) ) / h(i); % coarse mesh phi
    sigt(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 1 ));
    sigt(i)  = sigt(i)/(phi_c(i)*h(i)); % total cross-section
    sigs(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 5 ));    
    sigs(i)  = sigs(i)/(phi_c(i)*h(i));
    s_c(i)   = sum( dx(idx1:idx2)'.*s(idx1:idx2,1) ) / h(i);    
end
psi_c_ref(:,2) = Jp/mu_c;  % 2 is left to right!!
psi_c_ref(:,1) = Jm/mu_c;
phi_ref        = phi_c;

% COMPUTE d_sigs
NCM = length(phi_c);
% left-to-right
for i = 1:NCM
    exf = exp(-h(i)*sigt(i)/mu_c);
    tmp1 = 2*sigt(i)*( ...
           psi_c_ref(i+1,2) - psi_c_ref(i,2)*exf);
    tmp2 = 1 - exf;
    qbar = s_c(i) + sigs(i)*phi_c(i);
    if (fixed==1)
        d_sigs(i,2) =  ( tmp1/tmp2 - qbar ); 
    else
        d_sigs(i,2) = 1/phi_c(i) * ( tmp1/tmp2 - qbar );
    end
end
% right-to-left
for i = NCM:-1:1
    exf = exp(-h(i)*sigt(i)/mu_c);
    tmp1 = 2*sigt(i)*( ...
           psi_c_ref(i,1) - psi_c_ref(i+1,1)*exf);
    tmp2 = 1 - exf;
    qbar = s_c(i) + sigs(i)*phi_c(i);
    if (fixed==1)
        d_sigs(i,1) = ( tmp1/tmp2 - qbar );
    else
        d_sigs(i,1) = 1/phi_c(i) * ( tmp1/tmp2 - qbar );
    end
end


% here, we use an S2-like quadrature.  We want to conserve partial
% currents.  In S2, we have just one angle going right and one going left.
% The partial current J+ = w*mu*psi, or psi = J+/w/mu, and likewise for the
% leftward J-.  We correct Q such that psi(i+1)=J+/w/mu=f(psi(i),Q).

% We have reference partial currents going left and right
% convergence parameters
eps_phi = 1e-6; 
max_it  = 1000;
err_phi = 1;    
it = 0;
% Begin coarse mesh source iterations
ef = zeros(NCM,1);
for i = 1:NCM       % left-to-right
    ef(i)  = exp(-h(i)*sigt(i)/mu_c);
end
while (err_phi > eps_phi && it <= max_it )
    % Save old scalar flux
    phi0 = phi_c; 
    % Update sources
    for i = 1:NCM
        q_c(i) = s_c(i) + sigs(i)*phi_c(i);
    end
    % Perform sweeps (including the perturbed source)
    for i = 1:NCM       % LEFT-TO-RIGHT
        if (fixed==1) 
            tmp = d_sigs(i,2);
        else
            tmp = phi_c(i)*d_sigs(i,2);
        end
        psi_c(i+1,2) = 0.5*(q_c(i)+tmp)/sigt(i)*(1-ef(i)) ...
                       + psi_c(i,2)*ef(i);
    end
    for i = NCM:-1:1  	% RIGHT-TO-LEFT
        if (fixed==1)
            tmp = d_sigs(i,1);
        else
            tmp = phi_c(i)*d_sigs(i,1);
        end
        psi_c(i,1)   = 0.5*(q_c(i)+tmp)/sigt(i)*(1-ef(i)) ...
                       + psi_c(i+1,1)*ef(i); 
    end    
    % Update phi via quadrature & diamond difference relation (should use
    % analytic form...later...approximation should not affect basic idea)
    for i = 1:NCM
       phi_c(i) = 0.5.*(sum(psi_c(i,:))+sum(psi_c(i+1,:)));
    end
    % Update error and iteration counter
    err_phi =  max(  abs(phi_c-phi0)./phi0 );
    it = it + 1;
end

for i = 1:length(in.xfm)  
    % within cell indices
    idx1 = 1 + sum( in.xfm(1:(i-1)) ); % lower index
    idx2 = sum( in.xfm(1:(i  )) );     % upper index   
    phi(idx1:idx2,1)   = phi(idx1:idx2,1)  * phi_c(i)/phi_ref(i);
end
% hold on,figure(2),plot(phi_c./phi_ref,'k')
% lala=1;
end

function J = partcur( flag, psi, i, mu, w, in )
    if ( flag == 1 )
        idxs = in.ord/2+1:in.ord;
    else
        idxs = 1:in.ord/2;
    end
    J = sum( abs(mu(idxs)) .* w(idxs) .* psi(i,idxs)' );
end

